clear; clc;

StdVec      =   [1 1 1 -1 -1 60 -10]'*0.01/4;
PP          =   Setup_PP(struct('KAPPA',0.00/4,'XI',1/82.5/4,'Pir_SS',0.02/4,...
                                'Flag_dE',1,'Flag_Price',0,...
                                'B_Total_SS',0.91,'b_lb',-0.29,...
                                'TrProb_H',0.99,'Prob_H',0.37,...
                                'TrProb_ext',0.99,'Prob_ext',0.33,...
                                'Fin_DomShare',0.79,'Fin_AdjCost',0.8,...
                                'Cons_FracH',0.60,...
                                'RHO_M_dom',0.68,'RHO_M_ext',0.81,'RHO_Y_H',0.5,'RHO_p_F',0.90,...
                                'Taylor_Pi',1.1,'Taylor_ir',0.87,...
                                'Cons_ElasTN',6.19,'Cons_ElasHF',6.19));
%% Steady State

[SS,PP]     =   SteadyState_Solver(PP);

%% Construct the Model Structure
MODEL       =   Setup_MODEL(PP,SS);

% Check the Steady State in Equation Blocks
Res_NK      =   Equ_NK(PP,SS,0);
Res_Portfolio=  Equ_Portfolio(PP,SS,MODEL,0);
Res_HH      =   Equ_HH(PP,SS,MODEL,0);

if max(abs([Res_NK;Res_Portfolio;Res_HH]))<1e-5
    fprintf('Steady State is Correct!\n');
else
    error('Steaty State is not Correct, check it  ...');
end


%% Compute the Jacobian
EquJac          =   struct();
EquJac.HH       =   Equ_HH(PP,SS,MODEL,1,[],[],1);
EquJac.HH_Aux   =   Equ_HH_Aux(PP,SS,MODEL,1,[],[],1);
EquJac.NK       =   Equ_NK(PP,SS,1);
EquJac.Portfolio=   Equ_Portfolio(PP,SS,MODEL,1,[],[],1);

%% IRF

% Flexible Exchange Rate
PP.Taylor_dE    =   0;

EquJac.NK       =   Equ_NK(PP,SS,1);
fir_ord         =   JacAssemble(MODEL,EquJac);
tic;
[gx,hx,gxhx_ExitFlag,gxhx_Info]=gx_hx(fir_ord,1);
fprintf(['Perturbation is Solved, Elapse Time=',num2str(toc,'%.2f'),'s\n']);
if gxhx_ExitFlag~=1
    warning('Perturbation Solution is problematic');
    SOLUTION    =   [];
    return
else
    solution    =   struct('gx',gx,'hx',hx);
    
end
% StdVec          =   [1 1 1 -1 -1 60]'*0.01/4;
IRF             =   IRF_1order(MODEL,solution,StdVec,[],4*20);

SOLUTION        =   struct('PP',PP,'SS',SS,'MODEL',MODEL,'EquJac',EquJac,...
                           'fir_ord',fir_ord,'gxhx',solution,'IRF',IRF);
if ~strcmp(Options.FileName,'')
    save([Options.FileName,'.mat'],'PP','SS','MODEL','EquJac','fir_ord','solution','IRF');
end